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d 1 Introduction 

This paper considers a steady-state crack propagating along an interface between dissimilar 
anisotropic materials under an asymmetric load. As a model, a semi-infinite crack propagating 
at a constant speed along a perfect interface between two anisotropic half planes is analysed. 
Such a problem cannot be directly resolved using the approach developed by Ting (1980s) and 
Sih (1960s) by using Stroh formalism as was seen in Stroh (1950s). This approach should be 
adjusted to tackle non-symmetric loading. This has been done for the static problem in Morini 
(2012) where the Stroh formalism was combined with the weight function technique. In par- 
ticular, they found analytic expressions for the stress intensity factors. The main aim of this 
paper is to extend these results for a steady-state crack propagation and to propose a method 
for computing important fracture mechanics paramaters (e.g. stress intensity factors, energy 
\^ release rates) efficiently for an arbitrary loading that is non-smooth and does not necessarily 

00 have to be decreasing rapidly at infinity. We consider, as an example, only orthotropic materials 

while the approach can be extended to monoclinic materials as well. In addition, the second 
order asymptotic terms are also found for an arbitrary load, which is essential when performing 
perturbation analysis for the propagating interfacial crack. Finally, we present results of nu- 
^T) merical computations for a specific asymmetric load and discuss them in the context of fracture 

mechanics. 

> 

X 2 Preliminary Results 

In this section the mathematical model considered will be introduced along with the main 
mathematical results used for further analysis of the problem. A semi-infinite crack propagating 
at a constant speed, v, along a perfect interface between two semi-infinite anisotropic materials 
is considered. The crack is said to be occupying the region xi — vt < 0,X2 = as illustrated in 
Figure [T} 

The traction on the crack faces is known and is said to be given by the known functions 
o'2i{xi — vt,0^) = pf{xi — vt) for xi — vt < and body forces are assumed to be zero. No 
restrictions are imposed on the loading considered in this paper, that is: the loading need not 
vanish at ifinity or the crack tip. 

Section 2.1 introduces the previous results obtained for a static crack. These results include 
expressions for the traction along the interface and the displacement jump across the crack along 
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Figure 1: Geometry 



with their asymptotic expansions in the neighbourhood of the crack tip. Stress intensity factors 
are alos presented wlong with weight functions and the Betti formula for a static crack. Section 
2.2 discusses the dynamic crack and makes note of similar results as those presented for the 
static crack. Results that have been found for only the static problem that are going to be 
derived for the steady state in this paper are also noted. 



2.1 Static Crack 
2.1.1 Stroh Analysis 



The Stroh formalism, as used in Stroh (1962), is used in order to find expressions for both the 
displacement field and the stresses in the material. It makes use of Hooke's law, given by 



<7ij = Cijkieki = Cijki-g^, for i,j, k,l = 1,2, 



:i) 



where a is the stress, e is the strain and C is the stiffness tensor for the material. Furthermore, 
the following relationship relating the stress and displacement is also used 



E 

j,k,i=i 



dajj 
dxi 



0. 



Combining ([T]) and ^ it is seen that 



^ijkl^ Q — — U. 

OX j OX I 



(3) 



A solution for the displacement field is assumed in the form Ui = Aif(xi + px2) and this in 
conjunction with ([s]) gives the following eigenvalue problem 



(Cilfcl +pCilk2 +pCi2kl +p'^Ci2k2)Ak = 0. 
In matrix form this is written as 

(Q+p(R + R^)+p2T)A = 0, 



(4) 



(5) 
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where Q, R and T all depend on elastic constants for the given material. Solving this eigen- 
value problem to find the eigenvalues p and eigenvectors Ai the following relationships for the 
displacement vector and traction vector are obtained ( Suo , 1990 ) 



u(xi, X2) = Af (z) + Af (z) 



t(xi,X2) = Lf'(z) + Lf'(z), 



(6) 
(7) 



where z = xi + px2- Here the eigenvalues, p, are taken to be those with positive imaginary 
part and the matrix A is constructed from the eigenvectors of ([5]). The vector i{z) is given by 
ft = fi^i +PiX2)- The matrix L is given by 



L = (R^+pT)A. 

At this stage the matrix B = zAL~^ is introduced, where 
and W are also introduced, where 



-1. The bimaterial matrices H 



H = Bj + Bii, W = Bj - B,,, 

where the subscript / or II determines which material the matrix relates to. 

The analysis performed in Suo (1990) considered the same physical problem as shown in 
figure [T| with v = 0, with zero traction conditions imposed on the crack faces and continuous 
traction and displacement across the interface (xi > 0). The following results are found to be 
true along the real axis ( Suo[ 1990) 



h+(xi) +H"^Hh~(xi) = t(xi), < xi < 00, 



h+(xi) +H-iHh-(xi) 



0, 



-00 < xi < 0. 



(8) 



Here, h(z) is a function assumed in the form h.{z) = wz 2+**^. The branch cut of h{z) is placed 



along the negative real axis and both w and e are derived from the eigenvalue problem (Suo 

27re 



1990) 



Hw 



'Hw. 



An expression for the traction ahead of the crack was found in Suo ( 1990 ) 

1 



t(xi) 



^/2^^Xl 



Re{Kxfw), 



(9) 



(10) 



where w is found using the eigenvalue problem ([9]). K = Kj + iKjj is known as the stress 
intensity factor and includes the conributing factors to the traction of both mode / and mode 

II. 

The displacement jump across the crack, defined as [u], was also found in Suo (1990), where 
[.] is given by 

[/](x) = /(x,0+)-/(x,0-). (11) 



For xi < it was found that 
[u](xi) 



TT 



(H + H) 

cosh ire 



Re 



Ki-xiY'w 
l + 2ie 



(12) 



The asymptotic expansions of the physical traction and the jump in displacement across the 
interface as x — >• can be written as follows ( Morini , 2012[ ) 



[u](xi) = ^-^U{xi)K + ^-^^U{xi)Y2 + 4^W(xi)Y3 + 0((-xi)i: 
x/zTT \/2tt v27r 



(13) 
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t(xi) 



T{x,)K + ^r(xi)Y2 + -^r(xi)Y3 + 0{xl), 



(14) 



where K = [i^, -ftT] and Yj = [1^,1^]. Yi are constants derived in the same manner as the stress 
intensity factor K in order to find further terms in the asymptotic expansions. The matrices 
U{xi) and T{xi) are represented as follows 



Uixi) 



2(H + H) 

coshvre 



w 



{-xiY" w(-xi)- 



l + 2ie ' 1 - 2ie 



(15) 



2.1.2 Stress Intensity Factors 

An explicit formula for computing the stress intensity factor for symmetric loading was given in 



Suo (1990). It was shown that 



K 



2\2 /-o 

— coshvre / (—2:1) 2 *'^pi(xi)dxi, 
vr / ./_^ 



(16) 



where the vector pi(xi) is related to the applied traction p(xi) in the following way 



Pi 



w^Hp 
w^Hw 



2.1.3 Energy Release Rate 

A key component in the analysis of fracture mechanics is the determination of the energy release 
rate (ERR) when a unit area of interface is cracked. An expression was found for the ERR, 



denoted G, in Irwin (1957) 



](r)dr, 



(17) 



where A is an arbitrary length scale. It was found in Suo (1990) that using (10) and (12) that 



the following expression was found for the energy release rate 

w^(H + H)w|i^|2 



G 



4cosh^(7re) 



(18) 



2.1.4 Weight Functions 



The weight function U is now defined in the same vein as |WiUis| (|1995|). U = (C/i, U2f is the 



singular displacement field that is obtained in the problem where the static crack occupies the 
part of the axis with xi > is now considered. Therefore U is discontinuous over the positive 
portion of the real axis. The symmetric part of the weight function is given by [U](xi) and 
skew-symmetric part of the weight function is given by (U)(xi), where (.) is defined in the usual 
manner 



1 



{f){x) = -[f{x,^-^) + f{x,Q-)]. 



(19) 



The traction field associated with the displacement field, U, is denoted as Y = (Ti, ^2)^ and 
is said to be continuous over the interface (xi < 0) and the zero traction condition is imposed on 



the crack faces. Therefore, the following relationship stands for this problem, as seen in Morini 

h+(xi) + H~^Hh_(a;i) = for xi > 0, (20) 



(2012) 
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h+(xi) + H-^Hh_(xi) = T{xi] 
A solution for h{z) is assumed in the form 



for xi < 0. 



h(z) 



(21) 



(22) 



where the branch cut is now said to be along the positive xi-axis. This gives the eigenvalue 
problem 

(23) 



Hv = 



H is positive definite hermitian and therefore it is clear, by comparing (23) with i9h, that v = w 



. Therefore the singular traction in the steady state has the form (Morini, 2012) 



T(x) = ^^L^Re{R{-xif^), 
V 27r 



(24) 



where R = Ri + iR2 is a complex number in a similar fashion to the stress intensity factor for 
the physical problem. 

An expression relating the symmetric and skew-symmetric weight functions was found in 



Morini (2012) following from the work seen in Piccolroaz (2007) 

[U]^ix) = ^(isign(x)Im(H) -Re(H))t-(x), 



2\X 



■(isign(x)Im(W) - Re(W))T-(x), 



(25) 



(26) 



where the supercripts it denotes whether the function is analytic in the upper or lower half plane 
respectively. The " denotes the Fourier transform, defined as 



fix) 



f{x)e'^''dx. 



The Wiener-Hopf problem shown in (25) and (26) was solved in Morini (2012) and the 



corresponding parts of the weight function were shown there and will not be repeated in this 



paper for a static crack. However, later in this paper, equations (25) and (26) will be presented 



for the steady state despite the fact that they are obtained using the exact method used in 



Morini (2012). This is done because the representation used in this paper has a different branch 



cut to that in Morini (2012) which enables the calculation of further terms in the asymptotic 



expansions of the displacement and stress fields. Moreover, it is possible to obtain the equivalent 
result for a static crack by adopting the same method introduced for the steady state in this 
paper and setting v = 0. 



2.1.5 Betti Formula 

It was mentioned previously that there are now two displacement fields to consider; the physical 
displacment u and the singular displacement U. However, U is discontinuous for xi > whereas 
u is discontinuous for xi < 0. Also considered is the traction associated with U given by T 
which is continuous when xi < and the traction t associated with u which is continuous when 



xi > 0. It is possible to write an expression for the physical traction as (Morini, 2012) 



Cr(xi,X2)U2=0+ = p+(xi) +t(xi), O-(xi,X2)U2=0- = P- (xi) + t(xi). 



(27) 



Following the approach used in [Willis ( 1995| > and Piccolroaz (2007) the Betti formula is used 
to find a relationship between the weight functions and the physical solutions. It was seen that 



(28) 
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where 



7^ 




Here, (p) and [p] refer to the symmetric and skew-symmetric parts of the loading respectively. 



Further work peformed in Piccolroaz (2007) and Morini (2012), combining (25), (26) and 



(28), found an explicit expression for finding the stress intensity factor, K, using the weight 



functions and the loading applied on the crack faces. The following expression was obtained 



K 



27ri 



[U]+-(r)R(p)(T) + (U)^(r)R[p](r)dr, 



(29) 



where Z\ is a constant matrix derived from the asymptotic representation of (28). The matrix 



Z\ used in Morini (2012) uses different notation to that used in this paper but the same matrix 
can be found using results in section 3 with u = 0. 

In Suo (1990) symmetric loading was considered and therefore [p] = 0. It can easily be 



shown that both expressions for K, (16) and (29), are equivalent. See discussion in section 3. 



Following the method developed in Piccolroaz (2007) and Morini (2012) an expression for 



further asymptotic coefficients can be found depending on whether the applied loading is smooth 
and has a Fourier transform that vanishes at a fast enough rate at infinity. If this is the case 
the general expression for the asymptotic coefficients can be found using the equation 



1 



27ri 3 



/oo 
r^-H[U]+^(r)R(p)(r) + (U)^(r)R[p](r)}dT. 
-oo 



(30) 



Here, Zj is also derived from the asymptotic representation of (28) and is found for the steady 
state in section 3 of this paper and can be found for the static crack by setting t; = in the 
steady state formulation. 



2.2 Steady-State Case 

The steady state problem of figure[T]was considered in Yang (1991) and Yu (2000). In this paper 
it is assumed that the crack velocity v is constant and is lower than the Rayleigh wave speeds 
of both materials / and 11. 



2.2.1 Stroh Formalism and Suo Results for the Steady State 

When the crack is moving at the velocity v equation ^ becomes 

j,k,l=l ■' 



(31) 



where p is the material density. Using the same method as used for the static crack ([T]) and (31 ) 
are combined to give 

^'^''d^r''^- ^^^^ 

A new coordinate system is now introduced: (xi = xi — vt,X2 = X2). The following rela- 
tionship is therefore found in this new coordinate system 



ijkl 



d^Uk 
dxjdxi 



0, 



(33) 



where Cijki = Cijki - pv'^6ikdij6ii. 
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From this stage, for convenience, the steady state coordinates wih be written as xi = x and 
X2 = y- A solution is now assumed in the form Uj = Aif{x + py) the eigenvalue problem ^ is 
once again obtained. However, the matrix Q now relies on both the material constants and the 
crack velocity. The matrices R and T are identical to those in ([s]). 

Using the same methods as used for the static crack all equations from ([s]) to ( 16 ) still hold 
for the steady state with the coordinates {x, y). It is noted that the variables A, L, B, H, w and e 
depend on the crack velocity, as well as the elastic coefficents of the materials, for the steady 
state. 



2.2.2 Weight Function for the Steady-State Case 



The derivation of the weight function, U, for the steady-state coordinate system can be found 



using the same methods used for the static coordinate system, shown in Morini (2012). A differ 



ent representation is used in this paper allows a more convenient approach to the computation 
of higher order terms in the asymptotic expansions of the displacement and stress fields. 



Similarly to the static case seen in Morini (2012) the weight function, U, is obtained as the 



singular displacement field that is obtained in the problem where the crack occupies the part 
of the axis with x > is considered in order to obtain weight functions. The symmetric part 
and skew-symmetric part of the weight function are once again defined as [U](x) and (U)(x) 
respectively. The traction vector, denoted T, corresponding to this singular displacment field is 
assumed to be continuous over the interface (x < 0) and the traction free condition is imposed 



on the crack faces. Using the same methods as were used in Morini (2012) it can be shown that 



equations (25) and (26) still hold for the steady state. 



Expressions for the symmetric and skew-symmetric parts of the weight function, U, for 
orthotropic bimaterials are now found. To do this the associated traction, T, must be found. 



The method used for the static crack in Morini (2012) is used (with the results (36) and (37) 



written in a different way) . Substituting the value for w from ([TSj) into ( 24 ) , and saying R 
, yields the following result for — oo < x < 



1 



T^x) 




^[(-x)*^ + (-x)-^^] 



(34) 



Continuing to follow the method from Piccolroaz ( |2009 ) the same vectors is computed with 
R = i 

(-x)-i / -[(-x)- + (-x)-] 



THx) 




ll[(_3;)--(-x)- 



(35) 



The Fourier transforms of these two vectors is now computed. The branch cut for these 
vectors is situated along the positive real axis and polar coordinates with angle between —2tt 
and are taken. The identity 



r(fe + i) 

is also used. The Fourier transforms obtained are 



T^-(X) 



(ix)M / i[{-l-ie)T{l + ie){ixy 



for b > -1, 



ie)iixy 



(1 + 462) ^^gi [(-1 - ie)ril + ^e)iix^^ + (-1 + k)r(i - ze)(fx)^ 



(36) 




(l + 4e2)V^ 



ie){ix)-'' + {- 



ie)Ta+ie)iixy 



ieyixT'] 
- «e)(^x)* 



(37) 
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where r(.) is the gamma function and the branch cut of T is situated along the positive imag- 



inary axis. Note that the expressions (36) and (37) are written using a different representation 



than was used in Morini (2012). The reason behind this will become clearer in section 3. 



The Fourier transforms (25) and (26) can now be computed, for x G M, with the aid of 



equations (77) and (79). 



\x\ 
1 



-iPsign{xWHiiH22\ _ Hu 

^i/3sign(x)V^ii^22 J [ H22, 

-Hu 



-z/3sign(x) yjH 11H22 



\x\ U/3sign(x)\/-ffiii^22 



-H- 



T-(X), 



22 



(38) 



1 

1 

2W 



^-is\gT^{x)^/^^uH^n 
-i-fsign{x)VHiiH22 

-SiHii i7sign(x)\/i?iii?22\ 
-i7sign(x)\/-f/'ii-f^22 -(52-f^22 / 



6iHn 

, 52i^22, 



T-(X), 



(39) 



where branch cuts are now situated along the negative imaginary axis. 



2.2.3 Betti Formula for the Steady State 



It was shown in Willis (1995) that the Betti formula still holds for the steady state crack in 



isotropic materials. Using the same method it can be shown that the Betti formula still holds 
for the moving coordinate system in anisotropic materials. Therefore, the following expressions 
are found along the upper and lower parts of the real axis, respectively 



{U^(x' - X, 0+)7^c^(a;, 0+) - T^(x' - x, 0+)7^u(a;, 0+)}dx = 0, 



/oo 
{U^(x' - X, 0-)na{x, 0") - T'^(x' - X, 0")7^u(a;, 0")}dx = 0. 
-00 



(40) 



(41) 



The homogeneous case of ([s]) is now considered. Combined with the applied traction on the 
crack faces, p(x), the following expressions for traction are obtained 



o-2i(a;,y = 0+) = p+(x) + t(2;), C72i(a;,y = ) = p (x) + t(x). 



(42) 



Subtracting (41) from rt40| and using (27) along the definition of the symmetric and skew- 



symmetric parts of the weight function the following formula is obtained 

/oo 
{ [u]^(x' - x)nt{x) - r^{x' - X, o)7^[u] {x)}dx 
-00 

{[U]^(x' - x)n{p){x) + (U)^(x' - x)7^[p](a;)}dx, (43) 
Using the Fourier convolution theorem the following identity, which relates the Fourier trans- 



forms of the weight functions and the solutions of the physical problem, is obtained (Piccolroaz 
20071 ), ( |MoriEil[20T2l ) 



[u]+^m+ - f-^n[u]- = -[u]+^7^(p) - (u)^7^[p], (44) 

where the zt denotes whether the transform is analytic in the upper or lower half plane. 
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3 Evaluation of the Coefficients in the Asymptotic Expansion 
of the Displacement and Stress Fields 

3.1 Determination of the Stress Intensity Factor 

It is now possible to develop a method in order to find the stress intensity factor for orthotropic 
materials, similar to that seen for the static crack in Morini (2012). In the case of orthotropic 
materials the matrix T{x) in equation (14) is given by 

r(x) = 



-IX 



IX 



X 



(45) 



Note that this result is equivalent to (15) with the known value of w inserted. The Fourier 
transform of this expansion is computed in order to find the asymptotic expansion as x ~^ cxd, 
with Im(x) € (0, oo). The result is 



where 



-^^^i(x)K + 

-i(-ix)-"r(i + 



-ixY 
-ixY 



-^T(i + fe) 



-r(i 



''^^ '-T2(x)Y + 0((x)-i), 



i(-zx)"r(i-ie) ^ 
^(-ix)*T(i-ie) 



2V27r 



-ix)'T(j 



(46) 



(47) 



(48) 



It is noted here that these expressions differ to those seen in |Moririi| (2012) and Piccolroaz 



(2007) to incorporate the different branch cut used in this paper. It is now possible to find the 
asymptotic expansion of the members of Betti's identity from equation (44), using expressions 
( [SS] ) and ( 39 ) , as x — ^ oo 



= X^'^iK + x"'^2Y2 + x^'^sYg + 0{x-^), 
X-^ZiK + x"'22Y2 + x^'^sYg + 0(x"^), 



where Im(x) G (0,oo), (49) 
where Im(x) € (— oo,0). (50) 



The matrices Zi and Z2 are given by 

^11 



-l)(l-2ie) 
i(/3-l)(l-2t£) 



where 




^2(/3 + i)(i + 2ie) 
iE^(3 + l){l + 2ie) 



(/3-l)(l-2t£) (1+216) \ 

s+g- 1 

iE^il3+l){l+2ie) I ' 
s+g- / 



2r(^±ie)' 



2r(i±i6) 



Following the method of Morini ( 2012| ) ( |44| ) is written as a Riemann-Hilbert problem 

V'+(x) - V'-(x) = -[u]+^7^(p) - (u)^7^[p], (51) 

using the Plemelj formula it is possible to find ^l^{x) using the formula 



1 



^(r) 



dr. 



X 



(52) 
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The solution of this Riemann-Hilbert problem is given by 



where Im(x) S (0,oo), 
where Im(x) G (— oo,0). 



The asymptotic expansion of the Plemelj formula as x ~^ oo^ is given by 



1 

2TTi 



-2i 



,-3\ 



(53) 



Comparing the terms of this asymptotic expansion with the terms of the expansions (49) and 
(50) it is clear that = where Yi = K. Using (53) it is easily seen that the stress 

intensity factor, K, is given by 

1 r°° X 



K 



lim 

x-s>oo± 27ri 



[U]+^(r)7^(p)(r)-(U)^(r)7^[p](r) 



T-X 



dr, 



(54) 



where the explicit expression for Zi ^ is given by 



2s+s-(l + 4e^ 



£2 



(/3-l)(l-2i£) (/3_l)(l_2ie) 
1 i 



Assuming that the loading disappears in the region of the crack tip the limit in (54) exists and 
therefore the expression for the stress intensity factor, K, for the steady state agrees with that 
for the static case found in Morini (2012) (see equation (29)). 

Now that an expression for the stress intensity factro has been found it is possible to deter- 
mine an the energy release rate. As the velocities considered in this paper are all sub-Rayleigh 
they are all subsonic and therefore, as stated in Yu (2000), equation (17 can be used with 
the steady state-coordinates with A still an arbitrary length scale. Using (18) the following 
expression is obtained for the ERR in orthotropic materials 



G 



Hu{l-I3^)\K\ 



(55) 



3.2 General Expression for the Coefficients of the Higher Order Terms 



Using the asymptotic expansions ( 49 ) , ( 50 ) and the corresponding terms of and ( 53 ) a general 
expression for the jth coefficient of the asymptotic expansions, Yj, is found 



± 



lim 



X' 



2m{j - 1)! 



dJ 



-1 



X 



dr 



(56) 



This gives a general expression for the coefficents of the asymptotic expansion of the displacement 
and stress fields as 



lim 



± 2-Ki ■> 



r^■-l([U]+^(r)7^(p)(r) + (U)^(r)7^[p](r)) 



X 



X-T 



dr. 



(57) 



If the loading is applied in such a way that the limit exists it is easily seen that equation (57) is 
equivalent to (30). 



However, in general, the limit in (57) cannot be computed directly for j > 2. In this case an 
alternate method must be used in order to compute Yj. 
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Figure 2: Integration Shift in the x-Plane 



As the function p only exists the negative real x-axis its Fourier transform is analytic in the 
lower half x-plane. Therefore, [p] and (p) are also analytic in the lower-half plane. As long as 
the applied loading p vanishes within a region of the crack tip it is clearly seen that [p] and (p) 
decay exponentially as they tend to —ioo. It is also known that both [U]"'" and (U) are analytic 
in the lower-half plane apart from the negative imaginary axis. 

Using this infromation the contour of integration is moved to be either side of the negative 
imaginary axis, as shown in Figure [2j Equation (57) now becomes 



lim 



2-ni ^ 



X 



X 



dr 



X 



X ■ 



dr 



The limit of ( 58 ) can be taken to give 



1 



2TTi ^ 



ri-lf 



-V'(T)]dr, 



(5^ 



(59) 



where refers to the jump of the function i/j over the negative imaginary axis. 



The expression (59) can be simplified further by considering the continuity of (38) and (39). 



In order to do this ( 38 ) and ( 39 ) are rewritten in the following way 

prix) = 



{U){x) 




-i/3sign(x)\/i^ii-f^22 



i'ysign{x)V H11H22 

-i-fsign{x)VHuH22 



T-(X) 



T-(X)- 




The first term of both of these expressions simplify to give, respectively 





i(3^HnH22 



-i/3VHiiH22 




t (x), and 



1 







2x \i-fVfh^'. 



22 



-i^VHiiH22 




t-(x). 



Both of these terms are analytic in the lower half-plane and therefore continuous across the 
negative imaginary axis. For this reason they do not contribute to the general expression for 
the asymptotic coefficients, (59). 
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Figure 3: Loading 



A simpler expression for the asymptotic coefficients is now found 



27ri ^ 



where (j) is given by the expression 

(Pir) = Re(H)t-(r)7^(p)(r) +Re(W)t-(r)7^[p](r). 



(62) 



4 A Specific Example 

A specific example for computing the stress intensity factors for orthotropic materials is now 
considered. The loading on the crack faces is given by a point force of magnitude F acting on 
the upper crack face a distance a behind the crack tip and two point forces, both of magnitude 
F/2, acting on the lower crack face a distance b away from the point force acting upon the upper 
crack face. The loading moves at the same speed and in the same direction that the crack is 
propagating. This is shown in Figure |3j The forces are represented mathematically using the 
Dirac delta function (Piccolroaz 2009[ ) 



p+{x) = —F6{x + a), P-ix) 



F F 

--5{x + a + 6) - -5{x + a-b). 



(63) 



It is now possible to decompose the loading into its symmetric and skew-symmetric components 



{p){x) 



1 



[p+ix) +p-{x) 



IpKx) =P+{x) -p_(x) 



F F F 

-— 5(j; + a) — -^^{x + a — b) — -^^{x + a — 6), 

F F 

-FS{x + a) + —6{x + a + b) + —6{x + a-b). (64) 



In order to compute the stress intensity factors the Fourier transforms of the skew-symmetric 
and symmetric parts of the loading are required. These are given by 



{p)ix) 
\p\ix) 



F 



F 



-ix{a+b) 



F 



-ixia-b) 



__Pg-«Xa I _p-ix{a+b) _|_ _g-ix(a-6) 

2 2 



(65) 
(66) 
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It is now possible to compute expressions for the first and second order asymptotic coefficients, 



K and Y2, using expressions (54) and (57) respectively. 



To find an expression for K equation (54) is used. The solution is split into the parts 



corresponding to the symmetric and anti-symmetric parts of the loading, denoted and 
respectively 



H22 



(1 - /3) V Hii V TT 



F 



EH2 
(1-/3) 



H22 



-{1 + b/ay 



b/a) 



1 



(1-6/a) 



(67) 



(68) 



It is noted that for the loading considered in figure ([3]) that the limit in (57) does not exist 
and therefore equation (62) must be used with j = 2. Once again the coefficient is split into 



symmetric and anti-symmetric parts, given by 



Y4 



F 



(/3 - 1) V ^1 



H22 



\ + -Al + b/a) 



, E^52 
(/3-1) 



H22 



[I + b/a) 



-{I -b/a) 



--{1-b/ay 



(69) 



(70) 



Having computed expressions for the stress intensity factors it is now possible to calculate 
the energy release rate for two given materials. Material I is the piezoceramic Barium Titanate 
which has material properties Cn = 120.3GPa, C22 = 120.3GPa, Cu = 75.2GPa, Cqg = 21GPa 
and p = 6, 020kgm~^. Material II is set as Aluminium which has elastic constants Cn = 
107.3GPa,C22 = 107.3GPa,Ci2 = 60.9GPa,C66 = 28.3GPo and p = 2,700kgm-^. For the 
purpose of calculations, a is set as 1 in this paper. The velocity is normalised by dividing by 
cr, the lowest of the two Rayleigh wave speeds for the given materials, which for the materials 
considered is that of Barium Titanate. The enrgy release rate is non-dimensionalised in the 



results shown. G is normalised in the following manner: GC^l^/F'^. Here, C^g^ is taken as the 
value of Cee for the material above the crack. 

Fig. |4] shows the normalised ERR and Fig. [5] indicates the normalised energy release rates, 
and G^, corresponding to and K"^, respectively. Both components are normalised by 
the total energy release rate corresponding to K = + K^, given by G. 

The graph in Fig. |4] clearly shows that ERR increases as the velocity increases for all values 
of b/a and approaches infinity as the velocity approaches the Rayleigh wave speed , as expected. 
Moreover, the energy release rate is higher for larger values of 6/a, that is, as the asymmetry 
of the loading increases. Therefore, symmetric loading is energetically more beneficial than any 
choice of asymmetric one. 

Fig. [5] show that the contribution of the asymmetric part of the loading to the overall energy 
release rate is small compared to the contribution of the symmetric part for the lower crack 
velocities. However, when the crack is moving at velocities close to the Rayleigh wave speed 
the asymmetric nature of the loading begins to play a more crucial role in the calculations. 
This fact is also highlighted in the left hand graph of Fig. |6j Interestingly, the energy release 
rate corresponding to the symmetric part of the loading is higher than the overall ERR for 
the physical loading, apart from when the crack velocity is very close to the Rayleigh wave 
speed. The behaviour of the symmetric part of G near the Rayleigh wave speed is particularly 
noteworthy as it can be seen that for asymmetric loading (when b/a > 0) the ERR stops 
increasing and starts decreasing. Note that, for the symmetrical load {b/a = 0), the energy 
release does not depend on the crack speed. This was previously observed for isotropic and 



.(1) 



anisotropic bimaterials in (Wu, 1990 Yang 1991) 
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Figure 4: The normalised ERR, as a function of the velocity, for different positions of the 
self-balanced point forces applied to the crack surfaces, described by the ratio b/a in Fig. 1. 



Fig. shows the ratio of the mode 2 contribution of the stress intensity factor, K, to 
the mode 1 component. For the symmetric case the mode 2 component is for all values 
of the velocity, as one would expect, whereas for asymmetric symmetry the ratio is initially 
negative and then at a certain velocity the sign of the ratio changes. This is connected to 
the fact that for a determinate value of the propagation speed there is a change in the sign 
of the Dundurs parameter, /? (which only depends on the elastic constants of the materials). 
It is therefore possible to obtain a characteristic value for v, depending only on the material 
properties. This may lead to a possible change in the fracture mechanism or indicate a possible 
redirection of the interface crack propagation, for example kinking or branching. Note that 
there are other investigations, both theoretical and experimental, demonstrating an existence 
of a specific sub-Rayleigh velocity which is related to the stability of the crack propagation 
( |Obrezanova||2002b|aD . 
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Figure 5: a) The ratio /G, as a function of velocity, for different values ofb/a, with particular 
attention being given to the behaviour near c/j. b) The ratio G^/G with a different scale on the 
vertical axis. 
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Figure 6: a) The ratio G^/G^, function of the velocity, for different values of 6/a with 
particular attention being given to the behaviour near cr. b) The ratio K2/K1, as a function 
of thevelocity, for different values of b/a. c) The value of /3 as a function of the velocity. 
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A Orthotropic Stroh Matrices 

For orthotropic materials the matrices Q, R and T are given by 



Q 



Cii -pv^ 
. Cq&-pv\ 



,R 



Ci2 
^66 , 



,T 



Cm 
, C22, 



(71) 



Previously, expressions were found for the Stroh matrices for an orthotropic bimaterial with a 
crack propagating at speed v (Yang, 1991). The following paramaters are defined in the same 
manner as Yang (1991) 



a 



hp 



7/3 



Cm 



ai 



cTi 



02 



pv^ 
Cm 



and the variables 



iKii ^ On + '*ii'^22ai - (1 + ^^12) 

Oiia2\l and s = ^ 

K22 2aia2VKii'^22 
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It is seen that the eigenvalues, with positive imaginary part, of the equation 

(Q+p(R + R^)+/T)A = 



for the matrices from ( 71 ) are given by 



Pl,2 




(72) 



(73) 



Using the same normahsation as used in Yang (1991) the matrices A and L are given by 

A = 



1 

-Ai 



1 



Pi - Al l-p2A2^ 
Ki2 - K22P1A1 K22P2 - '^12A2 \ 



where 



A„ 



^ (1 + Ki2)Pf, 

It is now possible to find an expression for the hermitian matrix B given by 

B = ;AL~^ = ^ ( ^22Q2 ^2(1 + s)/C ^(^22 - Ki2ai/0\ 
CqgF y-i{K22 - Kual/O K22\/2^(l + s) J 

where F is the generalized Rayleigh wave function given by 

F = K22(k22? - 1 + a|) - nl^al/C 

It is possible to find the Rayleigh wave speed of a material by setting F = 0. 
The bimaterial matrix H can now be computed 

jj ^ / Hu -i(3VHuH22\ 
[il3^HnH22 H22 i 



From (76) it is seen that 
Hu = 



22 



At22aiv/2(l + g)/e 
CqqR 

' K22^j2i{l + S) 



CqqR 



+ 



+ 



At22QiV2(l + 5)/^ 
CqqR 

^122^2^(1 + 3) 



CqqR 



II 



H11H22 



K22 - K,VlG?2li 
CqqR 



II 



K22 - i^Yiajli 

CqqR 



(74) 
(75) 



(76) 



(77) 



In order to compute the weight functions the eigenvalues and eigenvectors of ([9]) are required. 
Using the representation (77) it is found that 



w 



1 [WL 

.2 V il22. 



27r 



In 



1-/? 
1 + /3 



(78) 
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Another key component for calculating the weight functions is the bimaterial matrix W 
Bj — B//. Using (76) it is seen that 



22 



(79) 



where 



K22-Kl2ai/5 


+ 


K22-Kl2ai/C 








II 



V-ff 11-^22 



2aiV2(l+s)/g 



2aiV2(l+s)/g ' 



77 



11 



ft^22V2g(l+s) 



2V2g(l+s) 



H, 



22 
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